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ABSTRACT 



We study the luminosity functions of high-redshift galaxies in detailed hydrodynamic simula- 
tions of cosmic reionization, which are designed to reproduce the evolution of the Lyman-a forest 
between z ~ 5 and z ~ 6. We find that the luminosity functions and total stellar mass densities 
are in agreement with observations when plausible assumptions about reddening at z ~ 6 are 
made. Our simulations support the conclusion that stars alone rcionizcd the universe. 

Subject headings: cosmology: theory - cosmology: large-scale structure of universe - galaxies: evolution 
- galaxies: formation - galaxies: high-redshift - galaxies: starburst - galaxies: stellar content - infrared: 
galaxies 



1. Introduction 

An important issue in cosmology today is the 
timing and mechanism of the reion ization of the 
universe. Since the early work of iMadau et al 



(|l999h . stars have been thought to be prime can- 
didates for the energy source, since the Gunn- 
Peterson trough disappears prior to the presumed 
epoch of quasars (z > 5). 

However, compounding the problem of iden- 
tifying the primary sources of cosmic reioniza- 
tion ar e recent HST observations of galaxies at 



z ~ 6 (IBouwens et all 120041: Bunker et al 



2004; 



Yan fc Wind horst! 12004 iGiayalisco et all 12004 : 
Hu et al.ll2005t IBouwens et al.ll2006h . It remains 
controversial whether enough stellar luminosity 
has actually been observed to account for reion- 
ization a t z ~ 6, as is seen in the spectra of SDSS 
quasars (White et al.ll2003l : iFan et alj|2006h . 

On t he othe r han d , both numerical simu- 
lations iGnedm] |2000t (Harford fc Gnedinl l2003t 



Springel k HernquistJ 120031: iGnedinl |2004 and 
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semi-analytical models (jNight et al. 200G) predict 
star formation rates at z ~ 6 nearly as large as 
those at z ~ 4, in dire ct conflict with som e of the 
observational studies ( Bunker et al.l 12004 ). but in 
agreement with recent st udies of Lyman -a emit- 
ting galaxies at z ~ 6 (|Hu et al. l2005lh which 
suggest that the star formation rate at z ~ 6 
is substantially higher than that found by HST 
observations. 

While the observational determination of the 
star formation rate at z ~ 6 will most likely re- 
main controversial for some time, in this paper we 
focus on a theoretical aspect of this problem. Of 
course, predicting a priori the absolute value of 
the average star formation rate in the universe at 
z ~ 6 is not possible with our current understand- 
ing of the galaxy formation process, in large part 
because the efficiency of star formation in molecu- 
lar clouds at high redshift is not known. However, 
modern cosmological simulations are sufficiently 
robust to predict the ratio of star formation rates 
at, say, z ~ 4 and z ~ 6. Thus, we can legitimately 
attempt to answer the following question: 

After imposing a requirement on a theoretical 
model to reproduce the observed evolution of the 
Lyman-a forest between z — 4 and z = 6, would 
it be possible, by appropriately adjusting the star 
formation efficiency as a free parameter, to fit si- 
multaneously galaxy luminosity functions at z w 4 
and z m 6 ? 
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In addition, we can use the recent observational 
estimates of the total stellar mass density at z w 5 
as an extra constraint. 

In this paper we use cosmological simulations of 
reionization as our theoretical model. The main 
advantage of our set of simulations is that they 
are not only able to resolve high redshift galaxies, 
but, by virtue of including radiative transfer of 
ionizing radiation, are also able to reproduce the 
mean properties of the obs erved Lyman-q fores t 
between redshifts 5 and 6 ( Gnedin fe Fan l2006h . 
In particular, our simulations are consistent with 
the end of the rei onization at z ~ 6, clearly seen 
in the SDSS data (|Fan et al]l200fih . 

Thus, we attempt to place both the observa- 
tions of the Lyman-a forest at z > 5 and the 
observations of star-forming galaxies at the same 
cosmic epoch into a unified picture, in which the 
same simulation is required to fit all available ob- 
servational data. 

We briefly describe the simulations in §2, 
present the star formation rates in §3.1, discuss 
numerical resolution issues in §3.2, present our 
luminosity functions in §3.3, and discuss the accu- 
mulation of stellar mass in §3.4. In §4 we conclude 
and discuss the results. 

2. Method 

Simulations used in this paper have been run 
with the " Softened Lagrangian Hydrodynamics" 
(SLH) code lGnedir] (j2000l 120041 ) . The simulations 
include dark matter, gas, star formation, chem- 
istry, and ionization balance in the primordial 
plasma, and three-dimensional radiative transfer. 
A flat ACDM cosmology with values of cosmo- 
logical param eters as determined by the first year 
WMAP data (jSpergel et al.ll2003f) is used through- 
out this paper. 

The radiative transfer is modeled using the Op- 
tically Thi n Variable Eddington T ensor (OTVET) 
method of iGnedin fe Abell pOOll ). While OTVET 
is an approximation, its main advantage is that it 
is self-consistently coupled to the rest of the sim- 
ulation code, and thus takes into account possible 
feedbacks among star formation, radiative trans- 
fer, and gas dynamics on spatial and temporal 
scales resolved in a simulation. 

Star formation is incorporated in the simula- 
tions using a phenomenological Schmidt law on 



scales below the resolution limit of the simulation 
(i.e. stars are allowed to form wherever the gas 
sinks below the resolution limit, irrespective of the 
large-scale environment or mass of the dark matter 
halo they are forming within) . This law introduces 
a free parameter: the st ar forma t ion effi ciency esF 
[as defined by eq. (1) of lGnedinl <|2000h ]. In order 
to incorporate the uncertainty in this parameter 
in our results, we use three different simulations 
in this paper. All of them have a simulation vol- 
ume of 8/i _1 comoving Mpc. The first two simu- 
lations include 128 3 dark matter particles and the 
same number of quasi-Lagrangian mesh cells for 
the gas evolution. The third simulation includes 8 
times more (256 3 ) resolution elements in the same 
simulation volume. 

There are two adjustable parameters in our sim- 
ulations: the efficiency of star formation esf and 
the ionizing efficiency, i.e. the amount of ionizing 
radiation emitted per unit mass of stars. This 
latter value depends upon the fraction of ioniz- 
ing photons escaping from the vicinity of a star 
forming region to the spatial scales resolved in the 
simulation, a value that is unknown and possibly 
variable. Our simulations make no assumptons 
about the value of the ionizing efficiencj0; rather 
the ionizing efficiency is adjusted to fit the ob- 
servational data on the Gunn-Peterson absorption 
in the spectra of S DSS quasars, as explained in 
Gnedin fc Far] (l2006h . 



The ionization history of the universe depends 
almost completely upon a product of the star 
formation efficiency and the ionizing efficiency, 
i.e. the factor that determines how many ioniz- 
ing photons are emitted per unit mass of dense 
and rapidly cooling gas (which is assumed to be 
transforming into stars). For each run the ion- 
izing efficiency is selected so that the simula- 
tion is consistent with the observed evolution of 
the mean tran smitted flux in the Lyman-a for- 
est at z - 6 fWhite et all 120031 : ISongaila|[2004 : 



Fan et al.ll2006T ). The star formation efficiency is 



then treated as the sole free parameter. 



1 Note that the ionizing efficiency parameter is related to 
the fraction of photons escaped from the resolution limit 
of a simulation. It is not directly related to the escape 
fraction from the virial radius, a quantity often used in 
analytical models. Computing the escape fraction from the 
virial radius is a formidable computational task, requiring 
much higher resolutuon simulations that we use here. 
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Table 1: Simulation Parameters 



Run 




AM b 






LowResHighSFR 


2 


2.6 x 10 Y 


0.15 


4 


LowResLowSFR 


2 


2.6 x 10 7 


0.05 


4 


HighResLowSFR 


0.64 


3.2 x 10 6 


0.05 


5.1 



"Spatial resolution in h~ 1 kpc. 
'Mass resolution in Mq. 
c Star formation efficiency. 
d Final redshift of simulation. 



The parameters of the three simulations are 
listed in Table [TJ The two smaller boxes (we label 
them "Low Res" ) are identical, except for the value 
of the star formation efficiency, which is different 
by a factor of 3 between the two simulations. The 
larger simulation ("HighRes") has the same star 
formation efficiency as the second LowRes simula- 
tion. To make the referral to a specific simulation 
clear, we use the names that reflect the resolution 
and star formation efficiencies of our simulations. 

Galaxies in the simulation are associated with 
gravitationally bound objects, identified with th e 
DENMAX algorithm (jBertschinger fc Gelbl[l99lh . 

Population synthesis is carried out using the 



STARBURST99 package (jLeitherer et al.l Il999h 
with the high m ass loss Geneva tracks as de- 
scribed previously ( Harford fc Gnedin 20031 ) using 
a metallicity of 0.001 (0.05 solar) and a galaxy es- 
cape fraction for ionizing photons of 0.1. The re- 
sults were adjusted to reflect a three exponent ini- 
ti al mass function, the "cano nical" function given 



Wcid ner fc Kroupal (|2006h . Also included is the 



emission of extra Lyman-a photons due to the 
recombination of ionized hydrogen. Throug hout 
the paper AB magnitudes ( Oke fc Gunn 19831 ) are 
used exclusively. 



Re ddening was computed according to lCalzetti 
(|l997l ). A color excess of 0.15 was included for 



z k, 4 galaxies ( Shaplev et al. 200ll ). For galaxies 
at z « 6 we present resu lts both without redden- 



ing (IStanwav et alj 120051) and with a color excess 
of 0.05 ( Bouwens et al.l I2006I ). For the luminos- 



ity functions in Figure [6] unreddened luminosities 
were multiplied by a factor of 1.8 to make them 
comparable to reddening with a color excess of 
0.05. 




Fig. 1 . — Star formation rate as a function of red- 
shift for our three simulations: HighResLowSFR 
(solid line), LowResLowSFR (dotted line), and 
LowResHighSFR (dashed line). The filled open 
circle shows the estimate of the averag e star for- 
mation rate at z = 4 from lSteidel et aL ( 19991 ). 



3. Results 

3.1. Star Formation Rate in the Simula- 
tions 

Figure [1] shows the evolution of the star for- 
mation history in the three simulations described 
above. As we emphasized above, all three runs 
are adjusted to reproduce the observed evolution 
of the mean transmitted flux in the Lyman-a tran- 
sition of neutral hydrogen above z ~ 5 reasonably 
well, so the reionization histories are very simi- 
lar in all three simulations and are consistent with 
the existing data. The star formation histories in 
the three runs are, however, quite different. The 
two low resolution simulations differ by approx- 
imately a factor of 3 by construction, since the 
star formation efficiency is different by that factor 
between them. The high resolution simulation has 
the same star formation efficiency as the LowRes- 
LowSFR run, but forms more stars at earlier times 
due to higher spatial resolution. In all simulations 
the star formation rate begins to flatten signifi- 
cantly at z < 7. 

We have examined in detail two snapshots in 
time: z — 4 (1.5 Gyr after the Big Bang), the 
latest time of the low resolution simulations, and 
z = 5.8 (0.99 Gyr after the Big Bang). Hereafter, 
unless specified otherwise, we use our HighRes- 
LowSFR simulation as a fiducial model for the 



z = 5.8 snapshot, and the LowResLowSFR sim- 
ulation as a fiducial model for the z = 4 snapshot; 
the HighResLowSFR simulation was not contin- 
ued until z — 4 due to computational expense. 
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Fig. 2. — Specific star formation rate (SFR per 
galaxy mass) vs the total galaxy mass for sim- 
ulated galaxies at z = 4 (top panel) and at 
z = 5.8 (bottom panel). The horizontal lines rep- 
resent galaxy mass ranges, and the vertical lines 
are standard deviations. The lines are omitted 
for the LowResLSFR case in the bottom panel 
for clarity. The dashed lin es are the results of 
Springel k Hernquist! (|2003l ). 



Deducing the global star formation rate by ob- 
serving individual galaxies requires knowledge of 
how star formation is distributed among galax- 
ies of different magnitudes. In order to test the 
robustness of the simulations in predicting this 
dependence theoretically, we compare our star 
formation rates in galaxi es of various masses 
with the simulations by ISpringel fc Hernquist! 
(|2003l ), who conducted extensive resolution and 



numerical convergence studies. The results of 
this comparison are shown in Figure [21 As one 
can see, the LowResLowSFR run underestimates 
specific star formation rates in galaxies with 
M < 2 x 10 10 M Q at both rcdshifts, while the 
HighResLowSFR is consiste nt with the results of 
Springel k Hernquist! (|2003h . 



It is important to underscor e that the physics of 
as co oling included in our and lSpringel k Hernquist 
2003) simulations is rather different: while 



Springel k Hernquist! ( 20031 ) only include equilib- 



rium cooling of primeval (i.e. metal-free) plasma, 
in our simulations non-equilibrium radiative trans- 
fer effects are taken into account. In addition, we 
include cooling from heavy elements and molecu- 
lar hydrogen. It is therefore somewhat unexpected 
that the two sets of simulations agree so well. On 
the other hand, if this agreement is not a mere 
coincidence, it emphasizes the robustness of mod- 
ern cosmological simulations in predicting total 
rates of star formation in well resolved galaxies 
(within the adopted phenomenological recipe of 
star formation, of course). 

It is interesting to note that for the same 
mass the higher redshift galaxies have a somewhat 
greater star formation rate than the lower redshift 
ones. Nevertheless the total star formation rate for 
these galaxies is greater at z = 4 than at z = 5.8 
because of the largest galaxies that have no coun- 
terpart at the earlier redshift. Since galaxies are 
observed only to limiting magnitudes, the mass 
evolution of galaxies by itself helps to explain why 
the luminosity function appears to e volve most at 
the bright end ( Bouwens et al. 20061 ). 

The star formation rates shown in Fig. [5] are 
measured in the simulations as the amount of stel- 
lar mass added in the previous 5 Myr. We have 
verified that, at z = 6, if that time interval is ex- 
tended up to 60 Myr, the results barely change. 
Thus, galaxies in the simulations have extended 
periods of star formation. This point is important 
because it suggests that the observed luminous UV 
galaxies which have been used to estimate the to- 
tal stellar density are a represent ative sample of all 



star-forming ga laxies at z ~ 6 (jStark et al. 2006; 



Yan et al J 1200a) . 
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3.2. Effects of Numerical Resolution 



Comparisons of apples and oranges are rarely 
useful, so in order to compare luminosities of sim- 
ulated and observed galaxies, we need to know 
whether observations count all the light emitted 
from individual resolved galaxies, or if a substan- 
tial fraction of light remains undetected below the 
sky surface brightness. If the former, then a com- 
parison between the simulations and the data is 
straightforward since it is easy to count all the 
emitted light in the simulations. 

If, however, a substantial fraction of the light is 
not detected in the observations, the comparison 
with the simulations becomes much more difficult. 
In this case a simulation is required not only to get 
the total luminosity of the galaxy right, but also to 
correctly model the full spatial distribution of this 
luminosity as well. H The latter presents a much 
more serious challenge to modern simulations. 

Unfortunately, we really have no way of know- 
ing how much light is missed in the observations. 
One possible way of approaching this problem 
could be by using sufficiently high resolution simu- 
lations to resolve star forming regions and recover 
the correct light profile. We, therefore, first need 
to discuss the issues of numerical resolution in our 
simulations, because in our simulations galaxies do 
hide a significant albeit not dominant fraction of 
their luminosity below the surface brightness limit 
of current observations. Thus, we need to under- 
stand whether this is merely a resolution effect, or 
a property real galaxies might have too. 

Numerical resolution affects simulated galaxies 
in diverse and complicated ways. Typically, simu- 
lations have constant spatial resolution in comov- 
ing units, with the real space physical resolution 
deteriorating with time. This effect would make 
later galaxies less concentrated than earlier ones. 
For our simulations, this effect is not tremendously 
large, since the difference between redshifts 4 and 
5.8 is only about 40%. In addition, the cool- 
ing consistency condition implemented in the SLH 



2 One may argue that to get the total luminosity right, a 
simulation has to resolve the star forming regions as well. 
This is not the case, however, since all star formation pre- 
scription used in modern simulations incorporate free ad- 
justable parameters, and one can phenomenologically fit 
the luminosity function of galaxies by properly adjusting 
parameters, without actually resolving the details of star 
formation. 
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Fig. 3. — Distributions of half light radii of sim- 
ulated galaxies brighter than zsso = 29.5 (unred- 
dened) at z y 6 co mpared to observed galaxies 
(jBunker et al.l I2004T ) . For our cosmology 1 arcsec 



is 5.94 proper kpc. 



code (|Gnedinl 119971 ) mitigates this effect further. 
Therefore, we might expect that if our spatial res- 
olution is sufficient at z = 5.8, then it is also suf- 
ficient at z — 4. 

However, one must also consider the mass res- 
olution, whose effect usually goes in the opposite 
direction. As time in the simulation progresses, 
galaxies become more massive on average, i.e. they 
are represented by a larger number of resolution 
elements (dark matter and stellar particles, and 
gas cells). The number of resolution elements 
in a given galaxy, if not sufficient, would affect 
the sharpness of the central peak, which in turn 
could affect the spatial distribution of star forma- 
tion. In our fiducial HighResLowSFR simulation a 
10 11 Mq halo consists of about 30,000 dark matter 
particles, a comparable number of baryonic cells, 
and about 10,000 stellar particles. Is this number 
sufficient to localize the star formation to within 
a few percent of the virial radius 1 ^ 

Figure [3] attempts to address this question by 
showing the distribution of half light radii for sim- 
ulated galaxies brighter than zgso = 29.5 (unred- 
dened). The high resolution galaxies are more 
compact than the low resolution ones and, in some 
cases, as compact as the observed ones, but there 



3 For reference, the comoving virial radius of a 10 11 NLq halo 
is about 100 comoving kpc, which translates into the an- 
gular diameter of 8" at z = 4 and 7" at z = 5.8. 
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Fig. 4. — Total reddened luminosity functions at 
z — 4 for the LowResLowSFR simulation (dashed 
line) based upon galaxies larger than 10 10 M Q . 
The solid black line with error bars shows the ob- 



servational data from|G^b^sch^t_ajJ |2004a|). The 
top panel shows the linear scale along the y-axis, 
while the bottom panel shows the log scale. The 
hatched band around the LowResLowSFR simu- 
lation shows Poisson errors. 

is still a substantial fraction (about 50%) of galax- 
ies that are larger than the observed ones. Also, 
the significant difference between the LowRes- 
LowSFR and HighRcsLowSFR runs indicates that 
we have not reached the numerical convergence on 
the half light radii of simulated galaxies. 

Therefore, in the rest of this paper, we will as- 
sume that our simulations do not properly resolve 
light profiles of most massive galaxies, and we will 
use the total luminosities of model galaxies. 

3.3. Luminosity Functions of Model Galax- 
ies 

It is customary to present the luminosity func- 
tion as a number density of galaxies per unit lumi- 
nosity, dn/dL, or the number of galaxies per unit 
magnitude, dn/dm. In this paper we, however, are 
primarily concerned with estimating the star for- 
mation rate at various redshifts from luminosity 
functions, so the quantity of interest to us is the 
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Fig. 5. — The same as in Fig. [5J but for z — 5.8, 
and showing in addition the HighResLowSFR run 
as the solid line within the hatched region. Two 
panels shows the unreddened and reddened lumi- 
nosity functio ns respectively. The observational 
data are from Bouwens et al. I (|2006l) . 



amount of light emitted per unit log in luminos- 
ity L 2 dn/dL, or, equivalently, per unit magnitude 
lCT - 4m dn/dm. 

Figures [4] and [5] present the main result of 
this paper. The total (i.e. including all the light 
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Fig. 6. — Comparison of faint ends of luminos- 
ity functions at z = 6. Shown are luminosity 
functions based upon the total galaxy popula- 
tion in the HighResLowSFR (red) and LowRes- 
LowSFR (green) simulations. Filled circles and 
triangles include all galaxies with total mass M > 
1.28 x 10 9 M (about 100 particles), and squares 
and diamonds all galaxies with total mass M > 
3.84 x 10 8 M Q (about 30 particles). Since the 
deepest observations extend only to the vertical 
dotted line, we have compared our results to the 
Schechter function fits to the observational data 



as presented in Table 12 of iBouwens et al.1 (12 0061 



2004h . 



The line sty les are long dash (IBouwens et ah 
20061). dots (iDickinson et alj l2004j) short dash 



(lYan & Win dhorstJl2004. dash dot feunker et al 



and dash dot dot (Malhotra et al 



20051 ). 



The luminosity functions are independent of the 
mass cut off out to at least zgso < — 15. 

of model galaxies irrespective of surface bright- 
ness) luminosity functions from all our simulations 
are compared with th e ob servational data from 
(|Gabasch et alJl2004al) and (jBouwens et al.ll2006h . 
We show luminosity functions both with linear and 
log vertical scale, because the log plot allows one to 
see the whole range spanned by the values, while 
the linear scale is useful because the total lumi- 
nosity (and, after an appropriate correction, the 
total star formation rate) is simply an area under 
the curve, which can be easily measured. 

The overlap in magnitudes between the simu- 
lations and the observations is not large; simula- 
tions, because of their limited box size, are not 
able to produce the brightest, most rare galaxies. 
Observations, on the other hand, are missing the 



Fig. 7. — The ratio of star formation rate to the 
stellar luminosity at 1500 A (as measured by the 
zsso magnitude) as a function of magnitude for 
the simulated galaxies. The luminosities are cal- 
culated from the unreddened magnitudes and plot- 
ted as a function of reddened magnitudes. Shown 
are galaxies larger than 10 10 M Q . The dashed line 
shows the generally assumed value for this ratio 
dMadau et al.lll99S 



majority of galaxies, that fall below the flux limit. 

The advantage of using our choice for repre- 
senting the galaxy luminosity function can now be 
illustrated. For example, it does appear that both 
the LowResLowSFR run and the observations are 
tracing the same luminosity function at z = 4. If 
we consider the combined curve as a true lumi- 
nosity function, it can be estimated from the top 
panel of Fig. 2] that observations (without correct- 
ing for incompleteness) account for about 50% of 
the total star formation rate, while the simulation 
includes about 60% of the total SFR (with the 
overlap between the simulation and observations 
being 10% of the total). 

The LowResHighSFR simulation (not shown) 
is a factor of 3 higher than the observational data, 
which is not surprising, since it has a 3 times 
higher star formation efficiency. Unless specifi- 
cally noted, we use the LowSFR simulations for 
comparisons to the observational data. 

The situation is different at z = 5.8. In the 
interval 27 < z§5o < 29, where the HighRes- 
LowSFR simulat ion overlaps with th e observa- 
tional data from Bouwens et al. (2006). the simu- 
lation predicts a factor of 2 more luminosity than 
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the data if no reddening c o rrecti on is included, as 
argued bv IStanwav et al. ( 2005f ). which is about 
a 3a deviation With the reddening correction of 
iBouwens et al.l (|2006l ). the simulation agrees with 
the data (accounting for about 70% of the total 
luminosity, while the data also account for about 
70% of the total light). The HighResLowSFR 
simulation is approximately consistent with the 
LowRcsLowSFR run, indicating that our simu- 
lated luminosity functions are close to the con- 
verged result (for the total luminosities of modeled 
galaxies). 

We, thus, underscore the crucial importance of 
knowing the reddening corrections at z ~ 6 with 
sufficient precision. Studies of reddening agree 
that it is significant l y less than at lower redshifts 



(sec 



Bouwens et al.l ([20061 ) for summary) 



A matter of considerable interest and debate is 
the slope of the faint end of the luminosity function 
when fitted to a Schechter function. Figure [^com- 
pares our luminosity functions to various fits to 
observational data from the literature. The verti- 
cal line indicates the faintest level observed, which 
i s ach ieved in the recent study bv IBouwens et al ' 



(120061) . We refer the reader to this paper (Fig. 15 
and accompanying discussion) for a comprehen- 
sive comparison of luminosity functions at z « 6 
from the literature. 

Simulation of the faint end of the luminosity 
function is limited by the minimum size galaxy 
that can be adequately resolved by the simula- 
tion. Since the exact number of dark matter 
and gas particles needed is uncertain, we show re- 
sults for lower cutoffs of 100 and 30 particles for 
the LowResLowSFR simulation. To meaningfully 
compare the HighResLowSFR run we use the cor- 
responding mass cutoffs (these will have 8 times as 
many particles). The results are the same out to 
z 850 < — 13. We find that, in contrast to the large 
galaxies of the previous figures, the smaller galax- 
ies tend to be considerably less efficient at star for- 
mation. For this reason, we think it unlikely that 
a simulation with greater mass resolution would 
change the faint end slope in the region brighter 
than Z850 < —15, fainter than can be currently 
observed. 

To convert from luminosity to star formation 
rate, it is usually assumed that a luminosity at 
1500 A of 8 x 10 27 erg/s/Hz correspo nds to 1 solar 
mass of newly formed stars per year (jMadau et al 
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Fig. 8. — The star formation history in simulated 
galaxies brighter than 29.5 magnitude corrected 
for the finite size of the simulation box. The up- 
per orange, hatched band show s the estimated 
range from Bouwens et al. ( 2006f) using the con - 



ventional conversion factor (Mada u et al. 1996). 



To facilitate comparison with the simulated galax- 
ies, which show a different conversion factor, the 
upper band has been transposed down by a factor 
of 1.5 (lower green band) as explained in the text. 



1998). This value, of course, depends upon the as- 
sumptions one uses for the population synthesis. 
We have used the mor e current "canonical" initia l 
mass function given in IWeidner fc Kroupal (120061 ). 
In addition, we use a metallicity of Z — 0.001 
(l/50th solar), which agrees better with the sim- 
ulation than does the usual solar metallicity as- 
sumption (most of the simulated galaxies have 
metallicities well below 10% solar.). Figure [7] 
shows that our conversion factors are about 1.5 
times higher than the conventional one. Of course, 
we do not know the true initial mass function and 
metallicity at these high redshifts. We emphasize 
the difference here to facilitate comparison of our 
results with the literature. 

Figure [51 which summarizes our findings, is 



based upon a similar figure from IBouwens et al 
(2006). The values are given to a limiting mag- 



nitude of 29.5. The upper hatched region shows 
the original range computed using the conven- 
tional conversion factor of 8 x 10 27 erg/s/Hz per 
lM Q /yr. To facilitate comparison with our re- 
sults, we also show this region transposed down- 
wards by a factor of 1.5 to form the lower hatched 
region, to correspond to the conversion factor of 
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1.2 x 10 28 erg/s/Hz per lM Q /yr from Fig. The 
star formation rates from our two LowSFR simula- 
tions with the same magnitude cut are also shown 
as individual symbols at z — 4 and z — 5.8. 
Both estimates agree with the data, although 
our value is s li ghtly lower at z — 4 than the 
Bouwens et al.l (12006 ) estimate. The small dis 



crepan cy is due to the fact that iBouwens et al 
(2006) estimate is based on the ISteidel et al 



1999) m easurement, while our simulations agree 



best with lGabasch et alj (|2004al ), who find a some- 
what lower val u e of st ar formation rate at z ~ 4. 
Gabasch et al. ( 2004bJ l considers their results to 



be in accord with those of ISteidel etHI (|l999l ). 



3.4. 



Total Accumulated Stellar Mass Den- 
sity 



Recently infrared observations have enabled es- 
timates of the accumulated stellar mass of high 
redshift galaxies. These results place constraints 
on prior generations of stars during the epoch of 
reionization and are an important test of cosmo- 
logical simulations. Much of the observational 
data involves galaxies too rare to be present in 
the simulation box. Even so, the trends in the 
data are interesting. 

Figure [9] shows a power law relation, which is 
nearly linear, between the star formation rate and 
the total stellar mass content of simulated galaxies 
at redshifts 4 and 5.8. These results agree broadly 
both in slope and magnitude with the hydrody- 
namic simulations of others ( Night et al.l 200' 



Springel & Hernquisd Eioi iFinlator et alJ I20O 



In addition, their data extends beyond our highest 
points to about 10 11 M Q . At redshif t 4 (top panel) 



the sm all overlap with the data of iFeulner et al 



(|2005l ) is encouraging because the trends are sim- 
ilar. The solid li ne showing good agre ement with 



the trend seen bv lFinlator et al. ( 20061 ) is notewor- 



thy because of the very different techniques used 
in their simulati ons. The semi-analytic model in 
lldzi et ah ( 20041 ). by contrast, has considerably 
more scatter at this red shift. The s pectro scopi- 
cally verified galaxies of Stark et al. (|2006l ) (bot- 



tom panel) exhibit a large spread, in contrast to 
the tight relationship in the simulated galaxies. 

Contrasting observational re sults at a r edshif t 
of about 6 have been found bv lYan et al. (2006). 
They find a population of galaxies with very high 
stellar masses but low star formation rates. These 
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Fig. 9. — Stellar mass as a function of the star 
formation rate at z — 4 (top panel) and z = 5.8 
(bottom panel) in the observations and the simula- 
tions (as labeled on the plots). Simulated galaxies 
shown are those larger than 1O 1O M . 

galaxies must have had higher rates of star forma- 
tion in the past since the current star formation 
rate is insufficient to have formed the estimated 
stellar mass in the time since the Big Bang. We do 
not see any counterpart to this population. How- 
ever, these galaxies are presumably very large and 
thus formed from density fluctuations too rare to 
be found i n our simu l ation box; the effective vol- 
ume of the Yan et al. ( 2006I ) observations is about 
600 times larger than the volume of our simulation 
box. 

The observational conversion from magnitude 
at 3.6 microns to stellar mass depends upon the 
unknown star formation histories of the galax- 
ies. Since we know the star formation history of 
the simulated galaxies, it is useful to compute di- 
rectly the expected observations. Fig. [TU] shows 



9 





0.5 




0.0 


c 




o 






-0.5 


£ 




to 






-1.0 


E 




1 

o 
in 


-1.5 


oo 




M 






-2.0 




-2.5 



A 

A 



25.5 26.0 



26.5 



27.0 
Z 850 



27.5 28.0 



28.5 



1.0 

ql 0.8 
u_ 
(/) 

o 

S 0.6 
u_ 

§ 0.4 

o 
o 

^ 0.2 

0.0 
0.0 



□ HIghResLowSFR 






J? A 


A LowResLowSFR 




ft 

- / ' A 


z = 5.803 : 



0.2 



0.4 0.6 
Time (Gyr) 



0.8 



1.0 




Fig. 10. — The (z 850 — 7713.6^) colors vs magnitudes 
in the zgso passband for the simulated galaxies at 
z = 4.0 (top panel) and z = 5.8 (bottom panel). 
The shaded band in the bottom panel shows the 
uppe r limits fo r the 79 IRAC invisible galaxies 
from lYan etafl (|2006l ). 



zs50 — 7n 3 . 6Al ) colors vs magnitude in the z850 pass- 
band for our two redshifts. There is little compa- 
rable observational data. We show in the bottom 



panel the IRAC invisible population of lYan et al 
1)20061 ) (Fig. Ell bottom panel). These are galax- 
ies, which, although detectable in the rest frame 
UV, have too little stellar mass to be detected in 
the optical (infrared at this redshift). Little can be 
said except that these limits are at least consistent 
with our simulation results. 

The star formation histories of the simulated 
galaxies larger than 10 10 M Q at z = 5.8 are remark- 
ably uniform. Figure [TT1 shows the average shape 
of this history to be a monotonically increasing 
function. Some individual galaxies show a decline 



Fig. 11. — The average star formation rate profile 
for simulated galaxies at z — 5.8 with total mass 
larger than 10 10 M0 as a function of time after 
the Big Bang. The error bars show the standard 
deviation. 

of star formation rate typically beginning at 0.8- 
1.0 Gyr. 

Figure [T2] shows the history of the total stellar 
mass density of the simulations. Note that the 
LowRcsHighSFR simulation accumulates about 
three times as much mass, as one would expect 
from the 3 times greater star formation efficiency. 
It is difficult to compare these results to obser- 
vations since these latter estimates are based on 
magnitude limited surveys. Nevertheless, these 
results reinforce our decision to use the lower star 
formation efficiency. The large arrow, s howing the 
estim ated stellar mass density at z = (jCole et al. 
2000), provides a ceiling, which should not be 
crossed. 

4. Conclusions and Discussion 

Despite its pivotal role in the evolution of the 
early universe, reionization remains poorl y under- 
stood . As has been previously argued (jGnedin 
2000), simulations provide important dynamical 
information inherently missing from analytic and 
semi-analytic approaches. This work is part of 
a continuing comparison of detailed cosmological 
simulations of reionization with observational data 
as they become available. 

We have shown that numerical simulations of 
reionization that reproduce the SDSS data on the 
evolution of Gunn-Peterson absorption in spectra 
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Fig. 12. — The evolution of the total stellar mass 
density in our three simulations (lines with sym- 
bols). The square, diamond, and square with X 
show observationa l resul ts of (iDrory et al. I l2005h . 
(jStark et alj l2006h . and (|Yan et aljbood) respec- 
tively. Als o shown for referen ce is the present day 
value from Cole et al. (j2000l) (arrow in upper left 
corner). 



of z ~ 6 quasars can also be made consistent with 
the observed galaxy luminosity functions at two 
different redshifts (z ~ 4 and z ~ 6) simultane- 
ously by adjusting a single parameter: a coefficient 
in the Schmidt law for star formation. The best-fit 
value of the coefficient is consistent with the obser- 



vational estimates from local galaxies (jKhokhlov 
1998). In the best-fit case, the star formation rates 
per unit galaxy mass are sim ilar to those found by 
Springel fc Hernquistj ( 2003 ) in their simulations 
using a very different numerical approach. 

The simulations and much of the analysis for 
this paper were completed be fore the third year 
WMAP results (|Sperge]|[2006h became available. 
The major change for structure formation is a low- 
ering of as- The expected result is a delay in 
star formation, but the exact effect, of course, is 
not clear without new simulations. We have em- 
phasized that the relative star formation rates at 
different redshifts are more amenable to analysis 
than are the absolute ones. Our single adjustable 
parameter would pro bably have to b e mod ified for 
the new cosmology. Alvarez et al. ( 20061) has es- 
timated the effect of the cosmological parameter 
changes on the collapsed fraction over a range of 
redshifts. The ratio of this fraction at z ~ 4 to 
that at z ~ 6 changes by a factor of about 1.5 



between the two cosmologies. The expected direc- 
tion of the effect would tend to bring the simu- 
lation into better agreement with observation (see 
Figure|5J). It should be noted that the other hydro- 
dynamic simulations to which we have compared 
our results have all used the higher value of a& . 

We are able to obtain agreement with the ob- 
served galaxy luminosity function at z ~ 6 if 
we consider all of the luminosity gravitationally 
bound to the galaxy and a dopt the reddening cor- 
rection of iBouwens et al' (2006). In addition, the 
total stellar mass densities at z « 5 to 6 are consis- 
tent with the latest observational estimates with 
the exceptions noted above. 

Our results suggest that Lyman Break galaxies 
observed at z sw 6 are the brightest ones at that 
time, with the caveat that the finite simulation 
box size limits the size of galaxies that we can 
simulate. 

The simulation argues that the amount of star 
formation we see is sufficient to reproduce the 
reionization history imprinted on the Lyman-a 
forest. 

The luminosity function of the simulation at 
z = 5.8 has been compared to the observational 
results of lBouwens et al. (2006) because their data 
overlap the simulation the most in magnitude. 
There seems to be general agreement that 28.5 
is the limiting magnitude at which completeness 
becom es an issue. Another group (jBunker et al 



2004) has taken the conservative approach of con- 
fining the analysis to this limit. They have argued 
for a lower star formation rate with the conse- 
quence that other reionization sources, in addition 
to stars, might be necessary. Our results are not 
easily reconciled with these estimates, even when 
the uncertainty in the free parameters in the sim- 
ulation, ic. the star formation efficiency and the 
effective escape fraction, are taken into account. 

iLe Fevre et al. ( 2005 ) have made an exhaus- 
tive study of all galaxies seen to a limiting magni- 
tude where the redshifts are determined spectro- 
scopically. They argue that a significant fraction 
of galaxies out to redshift 5 is missed by Lyman 
Break selection techniques and that these missed 
galaxies significantly affect the ionizing photon 
budget. Their magnitude limit is much too low 
for these galaxies to be seen in our simulations. 
However, if these conclusions are found to extend 
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to higher redshifts the simulations will have to 
be modified to get a realistic picture of reioniza- 
tion. Larger box sizes, which at present are pro- 
hibitively expensive, may be necessary. For the 
present paper it is important to note that their 
results do support the notion that stars are suffi- 
cient for reionization. 

All of our conclusions are limited by the size 
of the simulation box. The overall trend seen by 
Yan et al. I (l2006l) . in which the the largest galaxies 
have seen their zenith in the past, occurs over a 
much larger range of galaxy size than we can sim- 
ulate. Over smaller ranges there is no clear trend, 
and thus our results are not inconsistent with ex- 
isting data. However, we might have expected to 
see at least some trend in the simulation and w e 
do not. As discussed by iFinlator et al. ( 20061 ). 
the proportionality between star formation rate 
and stellar mass seen in hydrodynamic simulations 
may be an important limit ation of these mo dels. 
The semi- analytic model of lldzi et al. ( 2004 ). al- 
though it rests upon more assumptions, does at 
least produce more variability. Simulations with 
a much higher dynamic range may be needed to 
fully account for the population of galaxies con- 
tributing to reionization. 

Even if it turns out that other reionization 
sources, not included in our simulations, are neces- 
sary, the correspondence of the simulated galaxies 
to observations of this period is of considerable 
interest since the simulation appears at least to 
reproduce the ionization state of the universe be- 
tween redshifts 4 and 6. 
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